{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate raw combinations of elements or species \n",
    "This notebook aims to demonstrate the sheer scale of the combinatorial explosion. Based on some assumtions about how many elements (103) or species - elements in oxidation states - (403) there are, it calculates how many binary, ternary, quaternary... combinations there can be, within different integers stoichiometry limits.\n",
    "\n",
    "this notebook contains **ZERO** chemistry. \n",
    "\n",
    "The numbers that are produced match roughly to the raw combinations (before chemical filters) numbers obtained in Table 1 of the publication [Computational Screeningf of All Stoichiometric Inorganic Materials](https://www.cell.com/chem/fulltext/S2451-9294(16)30155-3) (D. W. Davies et. al, Chem, 2016).\n",
    "Small deviations from the numbers are due to changes in what we class as an \"accessible\" oxidation state for some elements.\n",
    "\n",
    "N.B. It can take a minute or two to calculate quaternaries on a desktop/laptop computer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports \n",
    "import itertools\n",
    "from fractions import gcd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In a search space of 403 elements:\n",
      "\n",
      "Number of unique combinations of 2 elements: 81003\n",
      "Number of unique combinations of 3 elements: 10827401\n",
      "Number of unique combinations of 4 elements: 1082740100\n",
      "\n",
      "Number of inequivalent binary ratios with max coefficient 2: \n",
      "3\n",
      "Number of inequivalent binary ratios with max coefficient 4: \n",
      "11\n",
      "Number of inequivalent binary ratios with max coefficient 6: \n",
      "23\n",
      "Number of inequivalent binary ratios with max coefficient 8: \n",
      "43\n",
      "Number of inequivalent ternary ratios with max coefficient 2: \n",
      "7\n",
      "Number of inequivalent ternary ratios with max coefficient 4: \n",
      "55\n",
      "Number of inequivalent ternary ratios with max coefficient 6: \n",
      "177\n",
      "Number of inequivalent ternary ratios with max coefficient 8: \n",
      "433\n",
      "Number of inequivalent quaternary ratios with max coefficient 2: \n",
      "15\n",
      "Number of inequivalent quaternary ratios with max coefficient 4: \n",
      "239\n",
      "Number of inequivalent quaternary ratios with max coefficient 6: \n",
      "1175\n",
      "Number of inequivalent quaternary ratios with max coefficient 8: \n",
      "3781\n",
      "\n",
      "Unique binary compounds (no stability constraints) with max coefficient 4: \n",
      "8.910e+05\n",
      "Unique binary compounds (no stability constraints) with max coefficient 6: \n",
      "1.863e+06\n",
      "Unique binary compounds (no stability constraints) with max coefficient 8: \n",
      "3.483e+06\n",
      "Unique ternary compounds (no stability constraints) with max coefficient 4: \n",
      "5.955e+08\n",
      "Unique ternary compounds (no stability constraints) with max coefficient 6: \n",
      "1.916e+09\n",
      "Unique ternary compounds (no stability constraints) with max coefficient 8: \n",
      "4.688e+09\n",
      "Unique quaternary compounds (no stability constraints) with max coefficient 4: \n",
      "2.588e+11\n",
      "Unique quaternary compounds (no stability constraints) with max coefficient 6: \n",
      "1.272e+12\n",
      "Unique quaternary compounds (no stability constraints) with max coefficient 8: \n",
      "4.094e+12\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dan/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:18: DeprecationWarning: fractions.gcd() is deprecated. Use math.gcd() instead.\n"
     ]
    }
   ],
   "source": [
    "# Functions\n",
    "def sum_iter(iterator):\n",
    "    return sum(1 for _ in iterator)\n",
    "\n",
    "\n",
    "def _gcd_recursive(*args):\n",
    "    \"\"\"\n",
    "    Get the greatest common denominator among any number of ints\n",
    "    \"\"\"\n",
    "    if len(args) == 2:\n",
    "        return gcd(*args)\n",
    "    else:\n",
    "        return gcd(args[0], _gcd_recursive(*args[1:]))\n",
    "\n",
    "\n",
    "def _is_irreducible(stoichs):\n",
    "    for i in range(1, len(stoichs)):\n",
    "        if gcd(stoichs[i-1], stoichs[i]) == 1:\n",
    "            return True\n",
    "    else:\n",
    "        return False\n",
    "\n",
    "# MAIN\n",
    "# >>>     Set to 103 for combinations of elements            <<<\n",
    "# >>> or 403 for combinations of element in oxidation states <<<\n",
    "search_space = 403\n",
    "\n",
    "\n",
    "def main():\n",
    "\n",
    "    compound_names = {2:'binary', 3: 'ternary',\n",
    "                      4: 'quaternary', 5:'quinternary'}\n",
    "\n",
    "    print (\"In a search space of {0} elements:\".format(search_space))\n",
    "    print (\"\")\n",
    "\n",
    "    def count_combinations(n):\n",
    "        c = sum_iter(itertools.combinations(range(search_space), n))\n",
    "        print (\"Number of unique combinations of {0} elements: {1}\".format(n, c))\n",
    "        return c\n",
    "\n",
    "    combinations = {n: count_combinations(n) for n in range(2,5)}\n",
    "\n",
    "\n",
    "    print(\"\")\n",
    "\n",
    "    # Count all the combinations of inequivalent ratios\n",
    "    # e.g. for ternaries up to x=2: (1,1,1), (1,1,2), (1,2,1),\n",
    "    #                               (1,2,2), (2,1,1), (2,1,2), (2,1,2)\n",
    "    def ineq_ratios_with_coeff(n_species, max_coeff):\n",
    "        print (\"Number of inequivalent {0} ratios with max coefficient {1}: \".format(\n",
    "                   compound_names[n_species], max_coeff)),\n",
    "        n_ratios = sum_iter(filter(\n",
    "            _is_irreducible,   # Filter operation cuts equivalent\n",
    "                               # stoichs e.g. accept (1,1), reject (2,2)\n",
    "            itertools.product(*( # All stoich combinations, e.g.\n",
    "                                  # Product of (1,2), (1,2) is\n",
    "                                  # (1,1), (2,1), (2,1), (2,2)\n",
    "                                  #\n",
    "                n_species * (tuple(range(1,max_coeff+1)),)\n",
    "                                    # Multiply nested tuples to form\n",
    "                                    # initial set of possibilities\n",
    "                                    # (1,2,..), (1,2,..),...\n",
    "                                ))))\n",
    "\n",
    "        print (\"{0}\".format(n_ratios))\n",
    "        return n_ratios\n",
    "\n",
    "    ratios = {n: {x: ineq_ratios_with_coeff(n, x) for x in range(2,9,2)} for n in range(2,5)}\n",
    "    print (\"\")\n",
    "\n",
    "\n",
    "    def unique_compounds(n_species, max_coeff):\n",
    "        print (\"Unique {0} compounds (no stability constraints) with max coefficient {1}: \".format(\n",
    "        compound_names[n_species], max_coeff)),\n",
    "        n_compounds = combinations[n_species] * ratios[n_species][max_coeff]\n",
    "        print (\"{0:5.3e}\".format(n_compounds))\n",
    "        return n_compounds\n",
    "\n",
    "    compounds = {n: {x: unique_compounds(n, x) for x in range(4,9,2)} for n in range(2,5)}\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    main()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
